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Abstract 

Second order buckling theory involves a one-way coupled coupled 
problem where the stress tensor from a plane stress problem appears 
in an eigenvalue problem for the fourth order Kirchhoff plate. In this 
paper we present an a posteriori error estimate for the critical buckling 
load and mode corresponding to the smallest eigenvalue and associated 
eigenvector. A particular feature of the analysis is that we take the 
effect of approximate computation of the stress tensor and also provide 
an error indicator for the plane stress problem. The Kirchhoff plate 
is discretized using a continuous/discontinuous finite element method 
which uses standard continuous piecewise polynomial finite element 
spaces which can also be used to solve the plane stress problem. 
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1 Introduction 


Buckling of thin plates can be modeled by an eigenvalue problem involving 
the stress tensor of the plane stress problem corresponding to a given load 
situation tangential to the plate. The smallest eigenvalues corresponds to 
the critical parameter multiplying the given plane stress load that results in 
buckling. 

Thin plates are modeled by fourth order differential equations accord¬ 
ing to the Kirchhoff-Love theory and require special attention when dis¬ 
cretized using the finite element method. In this paper we use the contin¬ 
uous/discontinuous Galerkin (c/dG) method proposed by Engel et al. [I] 
which is based on standard continuous piecewise polynomial spaces of order 
greater or equal to two inserted into a discontinuous Galerkin formulation, 
see Hansbo and Larson [3], of the fourth order plate equation. We refer 
also to Wells and Dung p] for a method closely related to the one presented 
here, and to Noels and Radovitzky [8] for an extension of the c/dG idea to 
Kirchhoff-Love shells. 

The c/dG formulation has the advantage that it uses standard finite el¬ 
ement spaces, is easy to implement, and extends naturally to higher order 
polynomials. Another important advantage in this particular problem is that 
we may solve the plane stress problem using the same finite element spaces. 
Note that this would not be the case if we, for instance, used nonconforming 
Morley elements for the plate problem since these element can not be used 
for the plane stress problem. 

In this paper we derive a posteriori error estimates for the critical buckling 
load and mode corresponding to the first eigenpair, and use these estimates 
to obtain mesh refinement strategies for error reduction. The error estimates 
are derived using duality techniques and are based on Larson [6] where a 
posteriori error estimates for the Poisson equation were presented. A partic¬ 
ular feature of the estimates presented herein is that we also take the effect 
of discretization of the plane stress problem into account. The error analysis 
of the buckling problem results in a specific goal functional which should be 
controlled in the plane stress solver. Here we follow the general approach 
to error estimation for one-way coupled problems developed by Larson and 
Bengzon [7], and adapted to linear second order plate theory in [2], We also 
mention the work [5] by Heuveline and Rannacher where a posteriori error 
estimates for a nonsymmetric eigenvalue problem related to the linearized 
stability of the Navier-Stokes equations is presented. These estimates also 
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involve the effect of the accuracy in the computed flow held on the eigenvalue 
problem and are thus related to our approach. 

This paper is organized as follows: in Section 2 we present the Kirchhoff- 
Love buckling problem and the continuous/discontinuous Galerkin method, 
in Section 3 we derive the a posteriori error estimates, in Section 4 we present 
some numerical results, and in Section 5 we present some conclusions. 


2 The Buckling Problem and Finite Element 
Method 

2.1 The Kirchhoff-Love Buckling Eigenvalue Problem 

The clamped Kirchhoff-Love buckling problem takes the form: find the plate 


displacements up (orthogonal to the plate) such that 

div div c P (V u p) — div t(cr M 'Vup ) = fp in 12 (1) 

Up = 0 on (912 (2) 

n ■ Ship = 0 on <912 (3) 

where t denotes the thickness of the plate and 

cr M = 2 pe(u M ) + Atr e{u M )I (4) 


where // and A are the Lame parameters, tr denotes the trace operator and 
I is the identity matrix, is determined by the membrane equation: find the 
membrane displacements % (tangential to the plate) such that 

-div <t m (u m ) = f M in ^ (5) 

Um = 0 on <912 (6) 

Here 

Et 3 

°~ P ( £ ) = 12(1 - u 2 ) (l 1 ~ U ) £ + ^(e) I) ^ 

is the plate stress tensor, e(v) = (Vv + (Vv) r )/2 is the strain tensor, E is the 
Young’s modulus, and v is the Poisson ratio, in terms of which A = E/(l + u) 
and /i = Ev/{ 1 — z/ 2 ). 
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Scaling the membrane forcing by a parameter A p, i.e., replacing the load 
by A pf M we note that by linearity cr M is replaced by Apcr M . The critical 
buckling loads are then determined by the eigenvalue problem: find Up and 
A p such that 


div div crp(Vup) — div t(Apcr M Vup) = 0 in Q (8) 

Up = 0 on dCl (9) 

n ■ Vup = 0 on dfl (10) 

The corresponding variational formulation reads: find the plate displace¬ 
ment Up G Hq(Q) and eigenvalue Ap G R such that 

ap(V«p, Vu) + \p(ct m X7up,Vv) = 0 Wv G Hq(Q) (11) 

where cr M G [L 2 (Cl)] 2x2 defined by flU, with u M G [//^(P)] 2 the solution of 

a M (u M ,v) = (f M ,v) Vv G [ffo(n)] J (12) 

Here the bilinear forms a P (-,-) and aw(-, ■) are defined by 

a P {0,-&) = (<T P (0),e(<?)) (13) 

a M {0, &) = {tr M (9), e(&)) (14) 

where (-, •) is the L 2 (Q) inner product. 

2.2 The Mesh and Finite Element Spaces 


We consider a subdivision T = {T} of into a geometrically conforming 
finite element mesh. We assume that the elements are shape regular, i.e., the 
quotient of the diameter of the smallest circumscribed sphere and the largest 
inscribed sphere is uniformly bounded. We denote by hp the diameter of 
element T and by h = maxpep hp the global mesh size parameter. We 
shall use continuous, piecewise polynomial, approximations of the transverse 
displacement: 


CVk = {ve C°(S2) : v\ T G V k (T) VT G T} (15) 

where Vk(T) is the space of polynomials of order k > 2 defined on T. Fur¬ 
thermore, we let CVk ,o = CVk n Hq. 
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We introduce the Scott-Zhang interpolation operator 7r : H ( ] (Q) —> C7\ 0 
and recall the following elementwise interpolation error estimate 

|U - TTu\ T ,m < Ch^lul^T^s (16) 

where 0<m<s<k + l and J\f(T) is the union of all elements which are 
neighbors to element T. 

To define our method we introduce the set of edges in the mesh, £ = {E}, 
and we split £ into two disjoint subsets 

£ = £j U £b (17) 

where £j is the set of edges in the interior of and £b is the set of edges 
on the boundary. Further, with each edge we associate a fixed unit normal 
n such that for edges on the boundary n is the exterior unit normal. We 
denote the jump of a function v G R at an edge E by [v] = v + — v~ for 
E G £i and [u] = v + for E e £ B , and the average ( v) = ( v + + v~)/2 for 
E G £i and ( v ) = v + for E € £b, where v ± = lim e | 0 v(x yen) with x e E. 

2.3 The Continuous/Discontinuous Galerkin Method 

We shall solve the membrane equation using standard continuous Galerkin 
and the plate problem with the continuous/discontinuous Galerkin method. 
The method takes the form: find Up G CVk P ,o an d Ap G R such that 

ApiyUp, Vv) + Apf(S M V[/p, Vv) — 0 Vv e CV kp ,o (18) 

where 'Em = 2f/,e(U m) + Atr e(U m)I with Um £ [CVk M ,o] 2 determined by 

a M {U M ,v ) = Vv G [CV kMfi f (19) 

The bilinear form Ap(-, •) is defined by 

A P {0,‘&)= J^(cr (0),e(0)) T 

tgt 

E(El£iL)£b 

- ([0], <ri - 

E^£j\J£b 

+ (2/i + 2A) 7 £ (20) 

Ee£j\J£b 
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for all 6, d e Q) TeT [H l (T)\ 2 . Here 7 is a positive parameter and hp is 
defined by 

h E = (\T + \ + \T~\)/(2\E\) for E = dT + n dT~ (21) 

with \T\ the area of T, on each edge E. See jlj for details on the value of 7. 

3 A Posteriori Error Estimates 

3.1 Preliminaries 

We first define a projector onto E(Xp) that is associated with the natu¬ 
ral scalar products involved in the variational statement. We define V\ p : 
H ( \ (Q) —» E(X P ) as follows 

A P (X7V\ p v,X7w) = A P (Vv,Vw) Vw e E(Xp) (22) 

Note that since w are eigenfunctions associated with Xp the projection also 
satisfies the following equation 

(<j m W\ p v, X7w) = (<7 M Vu, Ww) Vvj G E(Xp) (23) 

We introduce the norm 

IIMH 2 = A P (Vv, Vv), Vn 6 H 2 (Q) U CV kpfl (24) 

and normalize computed eigenvectors Up as follows 

IIIWpIH 2 = Ap(X/Up, X7U P ) = |Ap(E M Wp, W P )| = 1 (25) 

3.2 Error Representation Formulas 

The dual problem. To derive error representation formulas we introduce 
the following dual problem: find (j)p such that 

div div crp(V</)p) — div tAp((T A/ V0p) = Vfp bi H 

(pP = 0 on dVt 
n • V0p = 0 on dfl 
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(26) 

(27) 

(28) 


Different choices of the righthand side will lead to estimates for the errors in 
eigenvalues and eigenvectors. The righthand sides i/jp will be chosen in such 
a way that the solution to the dual problem is well defined. We return to 
these issues below. 

Multiplying with the error ep = Up — Up and integrating by parts we 
obtain 


(ep,ijjp) = (ep, div div crp(V0p)) 

- (e P , div tXp(cr mV0p)) 

= y ^(crp(Vep), £(V0 p))t 

rer 

- X] ([Ve P ],n • crp(V</>p))p 

E^SiUSb 

— Xpt(cr m^Up, Vcftp) 

— Ap(\7Up : V0p) — Apt(cr^Vf/p, V0p) (29) 

= — ApiyUp, v(0p — 7Tp0p)) 

— A pt(HM^Up, V(</>p — 7Tp(f>p)) 

+ (Ap — Xp)t(cr M VUp, V0p) 

— A pt((<r M — S M )Vt/p, V0p) (30) 

where in (l29]i we used the fact that [V0p] = 0 and then in (1301) we rearranged 
the terms using the identity Xpcr M = A P Yi M — (A P — Xp)cr M + Ap(<T M — 'E M ) 
and finally used Galerkin orthogonality (fTSll to subtract 7 T P (p P . 

Representation of the Error in the Eigenvalue Setting t/>p = 0 and 

denoting the solution to the dual problem by 0 p,a p we get 

(Ap — Xp)t(cr m^Up, V0p,Ap) 

= Ap(Vt/p, V (4>P,\ P — TTp^P.Ap)) 

+ Apf(£MVt/p, V(0p,Ap — np(j)p t \p)) 

+ A P t((Yi M — (Tm)'VUp, 'V(j)p t \ p ) (31) 

In this case the solution to the dual problem is an arbitrary eigenfunction 
associated with Ap, i.e. </>p,a p e E(X P ). Choosing </>p,a p = T , Apb r p/|||'PApb r p||| 
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we obtain the following estimate 


|(Ap — \ P )t(cr M 'VUp, V0p)| 

= |(Ap — Xp)Xp 1 A P (VUp, V0p iAp )| 

= |(A P - Xp)Xp l \A P (VV Xp U P , V0p, Ap )| 

= |(Ap-Ap)Ap 1 ||||ViP Ap 17p||| (32) 

We now assume that the computed eigenvalue Ap approximates the exact 
eigenvalue Ap and that there are constants 0 < 6 < 1 and h 0 such that 

|||V(/ — V\ p )Up\\\ < 6 (33) 

for all meshes with maxp e p hp < h 0 . We remark that the validity of this 
assumption follows from standard a priori convergence theory. Using (l32]h 
f[33|) . and the scaling (1251) together with Pythagoras identity we obtain 

|(A P — X P )t(cr M WUp, V0 P )| > |(Ap-Ap)Ap 1 |(l-<5 2 ) 1 / 2 (34) 

Finally, combining f[3T]) . (1341) . and using the triangle inequality we arrive at 

(1 — h 2 ) 1,/2 |(Ap — Ap)A P 1 | 

< \A P (yUp, V(</>p jAp — 7Tp0p, Ap )) 

+ Apt(S M VUp, V(0p, Ap — 7Tp0P, Ap ))| 

+ |Apt((<TM _ Sm)VUp, V0p, Ap )| (35) 

Representation of the Error in the Eigenvector. Following Larson 
[6] we define the error in an eigenvector to be the component orthogonal to 
the exact eigenspace which it approximates an element in. Note that this 
definition has the advantage that it covers also multiple eigenvectors. More 
precisely we will estimate the error in the H m (Q) seminorm for m = 0,1. We 
then define the error e m as 

e m = (I ~ Pm)Up (36) 

where P m is the orthogonal projection H m (Q) —> E(Xp) defined by (v — 
PoV,w ) = 0 and (V(i> — PiV, Vw) for all w G E(Xp) and m = 0,1, respec¬ 
tively. To represent the semi norm |ep| m we let V’p = V ’p,m = (—A) m ep iJn /|ep im 


with m — 0,1 and we denote the corresponding solution to the dual problem 
by 0p, m ,m = 0,1. We then get 

\ep,m\m = — Ap(VUp,'V((j)p ! m — KplfrP’m)) 

— Apt(SMVt/p, V(0P,n — TTp(j)p,m)) 

+ (Ap — Xp)t(cr M VUp , V0p, m ) 

+ Apt((SM — <Tm)Wp, V(j>p %m ) (37) 

In this case we require the solution 0p m to be orthogonal to E(X P ) to achieve 
uniqueness. 

Next we estimate the second term on the right hand side as follows 


Xpt(cr M 'VUp : 'V(j)p >rn ) — Ap(VU P , V0 P , m ) 

= Ap(Vep >m , V0p, m ) 
C | Cp,m | m 10P,m 14—m 
C Cm | Gp^m |m |0P,m | _m 
A C*m | Cp. rn | m 


(38) 


where we used the stability estimate 


10P,m |4—m A C'lV’P, m | —m 

and at last the identity | A Pm | -m = 1 which follows from the definition of 
-0 p, m ■ Thus we have 

|(Ap — Ap)f(crjvfVt/p, V0p, m )| 

< |(Ap — Ap)A p 1 |C m |ep im | m (39) 

Now again assuming that the computed eigenvalue Ap approximates the 
exact eigenvalue Ap and that there are constants 0 < S < 1 and h 0 such that 

\(A P -Xp)X- p 1 \C m <5 (40) 

for all meshes with maxp e p/rp < h 0 . We note again that the validity of this 
assumption follows from standard a priori convergence theory. Combining 
(|37jl . (J39J), and ()40l) and using the triangle inequality we obtain the estimate 

(1 — 5)|ep| m < |Ap(VCp, V(0p, m — 7Tp0p, m )) 

+ Apf(S M Vf/p, V(0 P , m — VTp0p im ))| 

+ |Apt((SM — (Tm)'VUp, V0p im )| (41) 
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Remark. The constant C m is of the form 


C m 


Cm 

gap(Ap) 


(42) 


where gap(Ap) is the distance between Ap and the closest eigenvalue. Thus 
assumption (140|) guarantees satisfactory resolution of the spectrum in the 
vicinity of Ap. 


Representation of the Modeling Error. Introducing the dual problem: 
find <Pm,x *= [tfoW su °h that 

a M (v, <f> M)X ) = Xt(cr M (v)VU P , V0p jJf ), (43) 

for all v G [Hq (fl)] 2 , X G {0,1, Ap}, we get, by setting v = e M := m m — 
Um an d using Galerkin orthogonality (ITTTj) for the membrane equation, the 
following error representation formula 

((Sm — 0 \fcf)VL/p, V(j)p t x) 

= (cr M (ejvf)Vt/p, V^ppc) 

= dM^M, 4 *M,x) 

= ClM( e iW) ^MM,X ) 

= (f Mi ^M^X ~ KM<t>M,x) ~ &m(UMi <Pm,X ~ n Al4 > M,X ) (44) 

3.3 Abstract A Posteriori Error Estimates 

Combining the estimates above we obtain the following abstract error esti¬ 
mates. For the error in the eigenvalue 

(1 — 5 2 ) 1 ^ 2 |(Ap — Ap)A P x | 

< \Ap(VUp, V(0p,A P — 7Tp0P,A P )) 

+ Apt(cr M -VL r p, V(0 p i a p — Tp<Ap,Ap))| 

+ I (/a/) 4>M,\p — ^A/^M.Ap) 

— CIm(UMi 0M,Ap — 7r A/</>M,Ap)l (45) 
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and for the error in the eigenvector 

(1 — S)\ep\ m < I Ap(VUp, V(0P,m — 7Tp0P,m.)) 

+ Apt(SMVf/p, V(0P, m — VTp0p jm ))| 

A l(/M) ^ 

— a M (U M , 0M,m — M,m) I (46) 

for m = 0,1. 

3.4 Error Estimates Using the Dual Weighted Resid¬ 
ual Approach 

Using standard procedures, involving integration by parts, the Cauchy-Schwartz 
inequality, a trace inequality, and the interpolation error estimate (fTFlh we 
obtain the following estimate 

\Ap(yUp , V (4>p t x ~ 7r P0p,x)) + 

A pt(TiM^Up,'V((j)p t x — 7Tp0p,v))| 

< ^ ^ Rp,t^Vp,x,t (47) 

TeTh 

where the plate element residual R PT and weight Wp t x,T are defined by 

Rp T — II/p — div div crp(WUp)\\% 

+ h^W [n • div crp(V[/p)]||| T 
+ ^ 3 ||[n- 0 -p(VUp)]||| T 
+ 7 2 hp 5 ||[VC/p]||| T (48) 

hUp,X,T = hrp P ’ X | (pp^X\j\f(K),ap i 0 < ttp^Y < kp + 1 (49) 

Here the regularity parameter a P) x reflects the regularity properties of the 
solutions to the dual problems. For the membrane problem we have the 
corresponding estimate 

\ a M(^KP 4 >m,x ~ 7T M < Pm,x)\ — Rm,tWm,X,T (50) 

TeT h 

where the residual and weight are defined by 

r 1i,t = II/m + div <tm(U m )\\t 

+ h^\\[n ■ div a M (U M )\\\ 2 aT (51) 
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and 


W M ,X,T — \<f>M,x\M{K),a M ,x ( 52 ) 

for 0 < oim,x < ku + G {Ap,0,1}. Collecting these estimates and the 
abstract a posteriori error estimates we finally arrive at the following dual 
weighted residual a posteriori error estimates 

(1 — <5 2 ) 1 / 2 |(Ap — Ap)A P 1 | < Rp,tWp,\,t 

T£T h 

+ Rm,tWm,\,t ( 53 ) 

TeT h 

and for the error in the eigenvector 

(1 — < 5 )|ep| m < Y^ Rp,tW p,m,T + Y^, RM,TWM,m,T ( 54 ) 

T&Th T&T h 

m — 0,1. Considering the expected optimal regularity of the dual problems 
we may expect 


«p,a — kp + 1 , ap.m — 4 — m, cum,™. — «m,a — 2 ( 55 ) 


3.5 Residual Based Estimates 

Using stability estimates for the solutions to the dual problems we obtain 
the residual based estimates 


(1 — <U)|(Ap — Ap)A P 1 
< 


C h T PXP R 2 p.T + Y h 

\TeT h TeT h 


2a.M,A p2 
T n M,T 


and for the error in the eigenvector 


(l-5) 2 |epl™< 


<c(j2 h£ p ' m R\ T + Y h r 

\T&T h T&T h 


' R 


M,T 


for m = 0,1. 


(56) 


(57) 
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4 Numerical examples 


4.1 Known stress tensor 

We consider the L shaped domain 

«:= (0,1) x (0,1)\ (1/2, l)x (0,1/2) 

The plate is simply supported on all boundaries (u — 0), and the in-plane 
stress tensor is chosen as the unit tensor. Thus, we have no error contribution 
from the membrane problem. We set E — 1, v — 1/4, and t — 1. We use 
the adaptive algorithm for the computation of the lowest three eigenvalues. 
The singularity in the inward-pointing corner is excited for the first two but 
not for the third, which is also clearly visible in the adaptation of the meshes 
shown in Figures 1-3. In Figure 4, we give the corresponding eigensolution, 
and in Figure 5 we give the corresponding effectivity indices (approximate 
error in eigenvalue divided by exact error). The third eigenvalue can be 
computed analytically, the first two have been estimated by an approximate 
solution on a dense mesh. The effectivity indices have been computed on a 
sequence of meshes obtained using a fixed ratio refinement technique where 
the elements with the highest 25% element error indicators have been refined 
in each step. The unknown constant in the error representation formula has 
been set so that the effectivity index is of medium size; the same constant 
has been used for all three eigenvalues. 

4.2 Computed stress tensor 

For our second example, we use the same domain, material data, and bound¬ 
ary conditions for the plate. For the elasticity computations, we use a body 
force f = (r, —9 r/10), where r denotes the distance from the inward point¬ 
ing corner. The boundary conditions were: clamped conditions at x — 1/2, 
y < 1/2, at y — 0, at x — 1, and at y — 1/2, x > 0. The remaining 
boundaries were traction free. 

In Figure 5 we give the adapted mesh using the full estimate, and, for 
comparison, we also give, in Figures 6-7, the corresponding meshes when 
only partial estimates, plate residual and stress residual, respectively, are 
used. In Figure 9 we show the lowest buckling mode for which the estimate 
is aiming. Finally, we show, in Figure 10, how the different residuals behave 
asymptotically as estimates of the eigenvalue error. Clearly, in order to 
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obtain an effectivity index that does not increase or decrease, we need the 
full residual, though we concede that the balance between these two residuals 
may be difficult to ascertain. We have here willfully chosen the balance in 
order to obtain a reasonably constant effectivity index for the full residual. 

5 Conclusions 

We have formulated a continuous/discontinuous Galerkin method for the thin 
plate buckling problem. The method has the advantage that we can solve 
both the membrane and plate problem with the same standard finite element 
spaces of continuous piecewise polynomials defined on triangles or quadri¬ 
laterals. Furthermore, we proved a posteriori error estimates for both the 
error in the eigenvalue (critical buckling load) and the eigenvectors (buckling 
modes) with the special feature that also the effect of approximate solution 
of the membrane problem is taken into account. Based on the estimates we 
constructed an adaptive algorithm for adaptive mesh refinement. 
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Figure 1: Adapted mesh for the first eigenvalue 
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Figure 2: Adapted mesh for the second eigenvalue 
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Figure 3: Adapted mesh for the third eigenvalue 
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Figure 5: Computed effectivity indices for the first three eigenvalue compu¬ 
tations. 




Figure 6: Adapted mesh using the full estimate. 
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Figure 7: Adapted mesh for a partial estimate (only the plate residual). 



Figure 8: Adapted mesh for a partial estimate (only the stress residual). 
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Figure 9: Buckling mode. 
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Effectivity index 



Refinement level 


Figure 10: Computed effectivity indices for the full estimate and the partial 
estimates (plate residual and stress residual). 
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